module cloud_rad_props

!------------------------------------------------------------------------------------------------
!------------------------------------------------------------------------------------------------

use shr_kind_mod,     only: r8 => shr_kind_r8
use ppgrid,           only: pcols, pver, pverp
use physics_types,    only: physics_state
use phys_buffer,      only: pbuf_size_max, pbuf_fld, pbuf_get_fld_idx, pbuf_old_tim_idx
use radconstants,     only: nswbands, nlwbands, idx_sw_diag, ot_length, idx_lw_diag
use abortutils,       only: endrun
use rad_constituents, only: iceopticsfile, liqopticsfile
use oldcloud,         only: oldcloud_lw, old_liq_get_rad_props_lw, old_ice_get_rad_props_lw, oldcloud_init

use ebert_curry,      only: scalefactor
use cam_logfile,      only: iulog

implicit none
private
save

public :: &
   cloud_rad_props_init,          &
   cloud_rad_props_get_sw,        & ! return SW optical props of total bulk aerosols
   cloud_rad_props_get_lw,        & ! return LW optical props of total bulk aerosols
   get_ice_optics_sw,             & ! return Mitchell SW ice radiative properties
   ice_cloud_get_rad_props_lw,    & ! Mitchell LW ice rad props
   get_liquid_optics_sw,          & ! return Conley SW rad props
   liquid_cloud_get_rad_props_lw, & ! return Conley LW rad props
   snow_cloud_get_rad_props_lw,   &
   get_snow_optics_sw

integer :: nmu, nlambda
real(r8), allocatable :: g_mu(:)           ! mu samples on grid
real(r8), allocatable :: g_lambda(:,:)     ! lambda scale samples on grid
real(r8), allocatable :: ext_sw_liq(:,:,:)
real(r8), allocatable :: ssa_sw_liq(:,:,:)
real(r8), allocatable :: asm_sw_liq(:,:,:)
real(r8), allocatable :: abs_lw_liq(:,:,:)

integer :: n_g_d
real(r8), allocatable :: g_d_eff(:)        ! radiative effective diameter samples on grid
real(r8), allocatable :: ext_sw_ice(:,:)
real(r8), allocatable :: ssa_sw_ice(:,:)
real(r8), allocatable :: asm_sw_ice(:,:)
real(r8), allocatable :: abs_lw_ice(:,:)

character(len=16) :: microp_scheme  ! microphysics scheme

real(r8) :: dmin = 2. *  13. / scalefactor
real(r8) :: dmax = 2. * 130. / scalefactor

! Minimum cloud amount (as a fraction of the grid-box area) to 
! distinguish from clear sky
! 
   real(r8) cldmin
   parameter (cldmin = 1.0e-80_r8)
!
! Decimal precision of cloud amount (0 -> preserve full resolution;
! 10^-n -> preserve n digits of cloud amount)
! 
   real(r8) cldeps
   parameter (cldeps = 0.0_r8)

! 
! indexes into pbuf for optical parameters of MG clouds
! 
   integer :: i_dei, i_mu, i_lambda, i_iciwp, i_iclwp, i_des, i_icswp

! indexes into constituents for old optics
   integer :: &
        ixcldice,           & ! cloud ice water index
        ixcldliq              ! cloud liquid water index


!==============================================================================
contains
!==============================================================================

subroutine cloud_rad_props_init()

   use netcdf
   use spmd_utils,     only: masterproc
   use ioFileMod,      only: getfil
   use error_messages, only: handle_ncerr
#if ( defined SPMD )
   use mpishorthand
#endif
   use phys_buffer,    only: pbuf_get_fld_idx
   use constituents,   only: cnst_get_ind
   use slingo,         only: slingo_rad_props_init
   use ebert_curry,    only: ec_rad_props_init, scalefactor

   character(len=256) :: liquidfile 
   character(len=256) :: icefile 
   character(len=256) :: locfn

   integer :: ncid, dimid, f_nlwbands, f_nswbands, ierr
   integer :: vdimids(NF90_MAX_VAR_DIMS), ndims, templen
   ! liquid clouds
   integer :: mudimid, lambdadimid
   integer :: mu_id, lambda_id, ext_sw_liq_id, ssa_sw_liq_id, asm_sw_liq_id, abs_lw_liq_id

   ! ice clouds
   integer :: d_dimid ! diameters
   integer :: d_id, ext_sw_ice_id, ssa_sw_ice_id, asm_sw_ice_id, abs_lw_ice_id

   liquidfile = liqopticsfile 
   icefile = iceopticsfile

   call slingo_rad_props_init
   call ec_rad_props_init
   call oldcloud_init

   i_dei    = pbuf_get_fld_idx('DEI')
   i_mu     = pbuf_get_fld_idx('MU')
   i_lambda = pbuf_get_fld_idx('LAMBDAC')
   i_iciwp  = pbuf_get_fld_idx('ICIWP')
   i_iclwp  = pbuf_get_fld_idx('ICLWP')
   i_des    = pbuf_get_fld_idx('DES')
   i_icswp  = pbuf_get_fld_idx('ICSWP')


   ! old optics
   call cnst_get_ind('CLDICE', ixcldice)
   call cnst_get_ind('CLDLIQ', ixcldliq)

   ! read liquid cloud optics
   if(masterproc) then
   call getfil( trim(liquidfile), locfn, 0)
   call handle_ncerr( nf90_open(locfn, NF90_NOWRITE, ncid), 'liquid optics file missing')
   write(iulog,*)' reading liquid cloud optics from file ',locfn

   call handle_ncerr(nf90_inq_dimid( ncid, 'lw_band', dimid), 'getting lw_band dim')
   call handle_ncerr(nf90_inquire_dimension( ncid, dimid, len=f_nlwbands), 'getting n lw bands')
   if (f_nlwbands /= nlwbands) call endrun('number of lw bands does not match')

   call handle_ncerr(nf90_inq_dimid( ncid, 'sw_band', dimid), 'getting sw_band_dim')
   call handle_ncerr(nf90_inquire_dimension( ncid, dimid, len=f_nswbands), 'getting n sw bands')
   if (f_nswbands /= nswbands) call endrun('number of sw bands does not match')

   call handle_ncerr(nf90_inq_dimid( ncid, 'mu', mudimid), 'getting mu dim')
   call handle_ncerr(nf90_inquire_dimension( ncid, mudimid, len=nmu), 'getting n mu samples')

   call handle_ncerr(nf90_inq_dimid( ncid, 'lambda_scale', lambdadimid), 'getting lambda dim')
   call handle_ncerr(nf90_inquire_dimension( ncid, lambdadimid, len=nlambda), 'getting n lambda samples')
   endif ! if (masterproc)

#if ( defined SPMD )
   call mpibcast(nmu, 1, mpiint, 0, mpicom, ierr)
   call mpibcast(nlambda, 1, mpiint, 0, mpicom, ierr)
#endif

   allocate(g_mu(nmu))
   allocate(g_lambda(nmu,nlambda))
   allocate(ext_sw_liq(nmu,nlambda,nswbands) )
   allocate(ssa_sw_liq(nmu,nlambda,nswbands))
   allocate(asm_sw_liq(nmu,nlambda,nswbands))
   allocate(abs_lw_liq(nmu,nlambda,nlwbands))

   if(masterproc) then
   call handle_ncerr( nf90_inq_varid(ncid, 'mu', mu_id),&
      'cloud optics mu get')
   call handle_ncerr( nf90_get_var(ncid, mu_id, g_mu),&
      'read cloud optics mu values')

   call handle_ncerr( nf90_inq_varid(ncid, 'lambda', lambda_id),&
      'cloud optics lambda get')
   call handle_ncerr( nf90_get_var(ncid, lambda_id, g_lambda),&
      'read cloud optics lambda values')

   call handle_ncerr( nf90_inq_varid(ncid, 'k_ext_sw', ext_sw_liq_id),&
      'cloud optics ext_sw_liq get')
   call handle_ncerr( nf90_get_var(ncid, ext_sw_liq_id, ext_sw_liq),&
      'read cloud optics ext_sw_liq values')

   call handle_ncerr( nf90_inq_varid(ncid, 'ssa_sw', ssa_sw_liq_id),&
      'cloud optics ssa_sw_liq get')
   call handle_ncerr( nf90_get_var(ncid, ssa_sw_liq_id, ssa_sw_liq),&
      'read cloud optics ssa_sw_liq values')

   call handle_ncerr( nf90_inq_varid(ncid, 'asm_sw', asm_sw_liq_id),&
      'cloud optics asm_sw_liq get')
   call handle_ncerr( nf90_get_var(ncid, asm_sw_liq_id, asm_sw_liq),&
      'read cloud optics asm_sw_liq values')

   call handle_ncerr( nf90_inq_varid(ncid, 'k_abs_lw', abs_lw_liq_id),&
      'cloud optics abs_lw_liq get')
   call handle_ncerr( nf90_get_var(ncid, abs_lw_liq_id, abs_lw_liq),&
      'read cloud optics abs_lw_liq values')

   call handle_ncerr( nf90_close(ncid), 'liquid optics file missing')
   endif ! if masterproc

#if ( defined SPMD )
    call mpibcast(g_mu, nmu, mpir8, 0, mpicom, ierr)
    call mpibcast(g_lambda, nmu*nlambda, mpir8, 0, mpicom, ierr)
    call mpibcast(ext_sw_liq, nmu*nlambda*nswbands, mpir8, 0, mpicom, ierr)
    call mpibcast(ssa_sw_liq, nmu*nlambda*nswbands, mpir8, 0, mpicom, ierr)
    call mpibcast(asm_sw_liq, nmu*nlambda*nswbands, mpir8, 0, mpicom, ierr)
    call mpibcast(abs_lw_liq, nmu*nlambda*nlwbands, mpir8, 0, mpicom, ierr)
#endif
   ! I forgot to convert kext from m^2/Volume to m^2/Kg
   ext_sw_liq = ext_sw_liq / 0.9970449e3_r8 
   abs_lw_liq = abs_lw_liq / 0.9970449e3_r8 

   ! read ice cloud optics
   if(masterproc) then
   call getfil( trim(icefile), locfn, 0)
   call handle_ncerr( nf90_open(locfn, NF90_NOWRITE, ncid), 'ice optics file missing')
   write(iulog,*)' reading ice cloud optics from file ',locfn

   call handle_ncerr(nf90_inq_dimid( ncid, 'lw_band', dimid), 'getting lw_band dim')
   call handle_ncerr(nf90_inquire_dimension( ncid, dimid, len=f_nlwbands), 'getting n lw bands')
   if (f_nlwbands /= nlwbands) call endrun('number of lw bands does not match')

   call handle_ncerr(nf90_inq_dimid( ncid, 'sw_band', dimid), 'getting sw_band_dim')
   call handle_ncerr(nf90_inquire_dimension( ncid, dimid, len=f_nswbands), 'getting n sw bands')
   if (f_nswbands /= nswbands) call endrun('number of sw bands does not match')

   call handle_ncerr(nf90_inq_dimid( ncid, 'd_eff', d_dimid), 'getting deff dim')
   call handle_ncerr(nf90_inquire_dimension( ncid, d_dimid, len=n_g_d), 'getting n deff samples')

   endif ! if (masterproc)

#if ( defined SPMD )
   call mpibcast(n_g_d, 1, mpiint, 0, mpicom, ierr)
!   call mpibcast(nswbands, 1, mpiint, 0, mpicom, ierr)
!   call mpibcast(nlwbands, 1, mpiint, 0, mpicom, ierr)
#endif

   allocate(g_d_eff(n_g_d))
   allocate(ext_sw_ice(n_g_d,nswbands))
   allocate(ssa_sw_ice(n_g_d,nswbands))
   allocate(asm_sw_ice(n_g_d,nswbands))
   allocate(abs_lw_ice(n_g_d,nlwbands))

   if(masterproc) then
   call handle_ncerr( nf90_inq_varid(ncid, 'd_eff', d_id),&
      'cloud optics deff get')
   call handle_ncerr( nf90_get_var(ncid, d_id, g_d_eff),&
      'read cloud optics deff values')

   call handle_ncerr( nf90_inq_varid(ncid, 'sw_ext', ext_sw_ice_id),&
      'cloud optics ext_sw_ice get')
   call handle_ncerr(nf90_inquire_variable ( ncid, ext_sw_ice_id, ndims=ndims, dimids=vdimids),&
       'checking dimensions of ext_sw_ice')
   call handle_ncerr(nf90_inquire_dimension( ncid, vdimids(1), len=templen),&
       'getting first dimension sw_ext')
   !write(iulog,*) 'expected length',n_g_d,'actual len',templen
   call handle_ncerr(nf90_inquire_dimension( ncid, vdimids(2), len=templen),&
       'getting first dimension sw_ext')
   !write(iulog,*) 'expected length',nswbands,'actual len',templen
   call handle_ncerr( nf90_get_var(ncid, ext_sw_ice_id, ext_sw_ice),&
      'read cloud optics ext_sw_ice values')

   call handle_ncerr( nf90_inq_varid(ncid, 'sw_ssa', ssa_sw_ice_id),&
      'cloud optics ssa_sw_ice get')
   call handle_ncerr( nf90_get_var(ncid, ssa_sw_ice_id, ssa_sw_ice),&
      'read cloud optics ssa_sw_ice values')

   call handle_ncerr( nf90_inq_varid(ncid, 'sw_asm', asm_sw_ice_id),&
      'cloud optics asm_sw_ice get')
   call handle_ncerr( nf90_get_var(ncid, asm_sw_ice_id, asm_sw_ice),&
      'read cloud optics asm_sw_ice values')

   call handle_ncerr( nf90_inq_varid(ncid, 'lw_abs', abs_lw_ice_id),&
      'cloud optics abs_lw_ice get')
   call handle_ncerr( nf90_get_var(ncid, abs_lw_ice_id, abs_lw_ice),&
      'read cloud optics abs_lw_ice values')

   call handle_ncerr( nf90_close(ncid), 'ice optics file missing')

   endif ! if masterproc
#if ( defined SPMD )
    call mpibcast(g_d_eff, n_g_d, mpir8, 0, mpicom, ierr)
    call mpibcast(ext_sw_ice, n_g_d*nswbands, mpir8, 0, mpicom, ierr)
    call mpibcast(ssa_sw_ice, n_g_d*nswbands, mpir8, 0, mpicom, ierr)
    call mpibcast(asm_sw_ice, n_g_d*nswbands, mpir8, 0, mpicom, ierr)
    call mpibcast(abs_lw_ice, n_g_d*nlwbands, mpir8, 0, mpicom, ierr)
#endif

    return

end subroutine cloud_rad_props_init

!==============================================================================

subroutine cloud_rad_props_get_sw(state, pbuf, &
                                  tau, tau_w, tau_w_g, tau_w_f,&
                                  diagnosticindex, oldliq, oldice)

! return totaled (across all species) layer tau, omega, g, f 
! for all spectral interval for aerosols affecting the climate

   ! Arguments
   type(physics_state), intent(in) :: state
   type(pbuf_fld),      intent(in) :: pbuf(pbuf_size_max)
   integer, optional,   intent(in) :: diagnosticindex      ! index (if present) to radiation diagnostic information

   real(r8), intent(out) :: tau    (nswbands,pcols,pver) ! aerosol extinction optical depth
   real(r8), intent(out) :: tau_w  (nswbands,pcols,pver) ! aerosol single scattering albedo * tau
   real(r8), intent(out) :: tau_w_g(nswbands,pcols,pver) ! aerosol assymetry parameter * tau * w
   real(r8), intent(out) :: tau_w_f(nswbands,pcols,pver) ! aerosol forward scattered fraction * tau * w

   logical, optional, intent(in) :: oldliq,oldice

   ! Local variables

   integer :: ncol
   integer :: lchnk
   integer :: k, i    ! lev and daycolumn indices
   integer :: iswband ! sw band indices

   ! optical props for each aerosol
   real(r8), pointer :: h_ext(:,:)
   real(r8), pointer :: h_ssa(:,:)
   real(r8), pointer :: h_asm(:,:)
   real(r8), pointer :: n_ext(:)
   real(r8), pointer :: n_ssa(:)
   real(r8), pointer :: n_asm(:)

   ! rad properties for liquid clouds
   real(r8) :: liq_tau    (nswbands,pcols,pver) ! aerosol extinction optical depth
   real(r8) :: liq_tau_w  (nswbands,pcols,pver) ! aerosol single scattering albedo * tau
   real(r8) :: liq_tau_w_g(nswbands,pcols,pver) ! aerosol assymetry parameter * tau * w
   real(r8) :: liq_tau_w_f(nswbands,pcols,pver) ! aerosol forward scattered fraction * tau * w

   ! rad properties for ice clouds
   real(r8) :: ice_tau    (nswbands,pcols,pver) ! aerosol extinction optical depth
   real(r8) :: ice_tau_w  (nswbands,pcols,pver) ! aerosol single scattering albedo * tau
   real(r8) :: ice_tau_w_g(nswbands,pcols,pver) ! aerosol assymetry parameter * tau * w
   real(r8) :: ice_tau_w_f(nswbands,pcols,pver) ! aerosol forward scattered fraction * tau * w

   !-----------------------------------------------------------------------------

   ncol  = state%ncol
   lchnk = state%lchnk

   ! initialize to conditions that would cause failure
   tau     (:,:,:) = -100._r8
   tau_w   (:,:,:) = -100._r8
   tau_w_g (:,:,:) = -100._r8
   tau_w_f (:,:,:) = -100._r8

   ! initialize layers to accumulate od's
   tau    (:,1:ncol,:) = 0._r8
   tau_w  (:,1:ncol,:) = 0._r8
   tau_w_g(:,1:ncol,:) = 0._r8
   tau_w_f(:,1:ncol,:) = 0._r8


   call get_liquid_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f)

   call get_ice_optics_sw   (state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, ice_tau_w_f)

   tau    (:,1:ncol,:) =  liq_tau    (:,1:ncol,:) + ice_tau    (:,1:ncol,:)
   tau_w  (:,1:ncol,:) =  liq_tau_w  (:,1:ncol,:) + ice_tau_w  (:,1:ncol,:)
   tau_w_g(:,1:ncol,:) =  liq_tau_w_g(:,1:ncol,:) + ice_tau_w_g(:,1:ncol,:)
   tau_w_f(:,1:ncol,:) =  liq_tau_w_f(:,1:ncol,:) + ice_tau_w_f(:,1:ncol,:)

end subroutine cloud_rad_props_get_sw
!==============================================================================

subroutine cloud_rad_props_get_lw(state, pbuf, cld_abs_od, diagnosticindex, oldliq, oldice, oldcloud)

! Purpose: Compute cloud longwave absorption optical depth
!    cloud_rad_props_get_lw() is called by radlw() 

   ! Arguments
   type(physics_state), intent(in)  :: state
   type(pbuf_fld),      intent(in)  :: pbuf(pbuf_size_max)
   real(r8),            intent(out) :: cld_abs_od(nlwbands,pcols,pver) ! [fraction] absorption optical depth, per layer
   integer, optional,   intent(in)  :: diagnosticindex
   logical, optional,   intent(in)  :: oldliq  ! use old liquid optics
   logical, optional,   intent(in)  :: oldice  ! use old ice optics
   logical, optional,   intent(in)  :: oldcloud  ! use old optics for both (b4b)

   ! Local variables

   integer :: bnd_idx     ! LW band index
   integer :: i           ! column index
   integer :: k           ! lev index
   integer :: ncol        ! number of columns
   integer :: lchnk

   ! rad properties for liquid clouds
   real(r8) :: liq_tau_abs_od(nlwbands,pcols,pver) ! liquid cloud absorption optical depth

   ! rad properties for ice clouds
   real(r8) :: ice_tau_abs_od(nlwbands,pcols,pver) ! ice cloud absorption optical depth

   !-----------------------------------------------------------------------------

   ncol = state%ncol
   lchnk = state%lchnk

   ! compute optical depths cld_absod 
   cld_abs_od = 0._r8

   if(present(oldcloud))then
      if(oldcloud) then
         ! make diagnostic calls to these first to output ice and liq OD's
         !call old_liq_get_rad_props_lw(state, pbuf, liq_tau_abs_od, oldliqwp=.false.)
         !call old_ice_get_rad_props_lw(state, pbuf, ice_tau_abs_od, oldicewp=.false.)
         ! This affects climate (cld_abs_od)
         call oldcloud_lw(state,pbuf,cld_abs_od,oldwp=.false.)
         return
      endif
   endif

   if(present(oldliq))then
      if(oldliq) then
         call old_liq_get_rad_props_lw(state, pbuf, liq_tau_abs_od, oldliqwp=.false.)
      else
         call liquid_cloud_get_rad_props_lw(state, pbuf, liq_tau_abs_od)
      endif
   else
      call liquid_cloud_get_rad_props_lw(state, pbuf, liq_tau_abs_od)
   endif

   if(present(oldice))then
      if(oldice) then
         call old_ice_get_rad_props_lw(state, pbuf, ice_tau_abs_od, oldicewp=.false.)
      else
         call ice_cloud_get_rad_props_lw(state, pbuf, ice_tau_abs_od)
      endif
   else
      call ice_cloud_get_rad_props_lw(state, pbuf, ice_tau_abs_od)
   endif
      
   cld_abs_od(:,1:ncol,:) = liq_tau_abs_od(:,1:ncol,:) + ice_tau_abs_od(:,1:ncol,:) 

end subroutine cloud_rad_props_get_lw

!==============================================================================

subroutine get_snow_optics_sw   (state, pbuf, tau, tau_w, tau_w_g, tau_w_f)
   type(physics_state), intent(in) :: state
   type(pbuf_fld),      intent(in) :: pbuf(pbuf_size_max)

   real(r8),intent(out) :: tau    (nswbands,pcols,pver) ! extinction optical depth
   real(r8),intent(out) :: tau_w  (nswbands,pcols,pver) ! single scattering albedo * tau
   real(r8),intent(out) :: tau_w_g(nswbands,pcols,pver) ! assymetry parameter * tau * w
   real(r8),intent(out) :: tau_w_f(nswbands,pcols,pver) ! forward scattered fraction * tau * w

   real(r8), pointer, dimension(:,:) :: dei, iciwpth
   real(r8) :: dlimited ! d limited to dmin,dmax range

   integer :: i,k,i_d_grid,k_d_eff,i_swband,ncol,lchnk
   real :: wd, onemwd, ext, ssa, asm

   ncol = state%ncol
   lchnk = state%lchnk

! temporary code to support diagnostics of snow radiation
   iciwpth => pbuf(i_icswp)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
   dei     => pbuf(i_des)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
! temporary code to support diagnostics of snow radiation

   do i = 1,ncol
      do k = 1,pver
         dlimited = dei(i,k) ! min(dmax,max(dei(i,k),dmin))
         if( iciwpth(i,k) < 1.e-80_r8 .or. dlimited .eq. 0._r8) then
         ! if ice water path is too small, OD := 0
            tau    (:,i,k) = 0._r8
            tau_w  (:,i,k) = 0._r8
            tau_w_g(:,i,k) = 0._r8
            tau_w_f(:,i,k) = 0._r8
         else 
            !if (dlimited < g_d_eff(1) .or. dlimited > g_d_eff(n_g_d)) then
               !write(iulog,*) 'dei from prognostic cldwat2m',dei(i,k)
               !write(iulog,*) 'grid values of deff ice from optics file',g_d_eff
               !call endrun ('deff of ice exceeds limits')
            !endif
            ! for each cell interpolate to find weights and indices in g_d_eff grid.
            if (dlimited <= g_d_eff(1)) then
               k_d_eff = 2
               wd = 1._r8
               onemwd = 0._r8
            elseif (dlimited >= g_d_eff(n_g_d)) then
               k_d_eff = n_g_d
               wd = 0._r8
               onemwd = 1._r8 
            else
               do i_d_grid = 1, n_g_d
                  k_d_eff = i_d_grid
                  if(g_d_eff(i_d_grid) > dlimited) exit
               enddo
               wd = (g_d_eff(k_d_eff) - dlimited)/(g_d_eff(k_d_eff) - g_d_eff(k_d_eff-1))
               onemwd = 1._r8 - wd
            endif
            ! interpolate into grid and extract radiative properties
            do i_swband = 1, nswbands
               ext = wd*ext_sw_ice(k_d_eff-1,i_swband) + &
                 onemwd*ext_sw_ice(k_d_eff  ,i_swband) 
               ssa = wd*ssa_sw_ice(k_d_eff-1,i_swband) + &
                 onemwd*ssa_sw_ice(k_d_eff  ,i_swband) 
               asm = wd*asm_sw_ice(k_d_eff-1,i_swband) + &
                 onemwd*asm_sw_ice(k_d_eff  ,i_swband) 
               tau    (i_swband,i,k)=iciwpth(i,k) * ext
               tau_w  (i_swband,i,k)=iciwpth(i,k) * ext * ssa
               tau_w_g(i_swband,i,k)=iciwpth(i,k) * ext * ssa * asm
               tau_w_f(i_swband,i,k)=iciwpth(i,k) * ext * ssa * asm * asm
            enddo
         endif
      enddo
   enddo

   return
end subroutine get_snow_optics_sw   

!==============================================================================
! Private methods
!==============================================================================

subroutine get_ice_optics_sw   (state, pbuf, tau, tau_w, tau_w_g, tau_w_f)
   type(physics_state), intent(in) :: state
   type(pbuf_fld),      intent(in) :: pbuf(pbuf_size_max)

   real(r8),intent(out) :: tau    (nswbands,pcols,pver) ! extinction optical depth
   real(r8),intent(out) :: tau_w  (nswbands,pcols,pver) ! single scattering albedo * tau
   real(r8),intent(out) :: tau_w_g(nswbands,pcols,pver) ! assymetry parameter * tau * w
   real(r8),intent(out) :: tau_w_f(nswbands,pcols,pver) ! forward scattered fraction * tau * w

   real(r8), pointer, dimension(:,:) :: dei, iciwpth
   real(r8) :: dlimited ! d limited to dmin,dmax range

   integer :: i,k,i_d_grid,k_d_eff,i_swband,ncol,lchnk
   real :: wd, onemwd, ext, ssa, asm

   ncol = state%ncol
   lchnk = state%lchnk

   iciwpth => pbuf(i_iciwp)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
   dei     => pbuf(i_dei)%fld_ptr(1,1:pcols,1:pver,lchnk,1)

   do i = 1,ncol
      do k = 1,pver
         dlimited = dei(i,k) ! min(dmax,max(dei(i,k),dmin))
         if( iciwpth(i,k) < 1.e-80_r8 .or. dlimited .eq. 0._r8) then
         ! if ice water path is too small, OD := 0
            tau    (:,i,k) = 0._r8
            tau_w  (:,i,k) = 0._r8
            tau_w_g(:,i,k) = 0._r8
            tau_w_f(:,i,k) = 0._r8
         else 
            !if (dlimited < g_d_eff(1) .or. dlimited > g_d_eff(n_g_d)) then
               !write(iulog,*) 'dei from prognostic cldwat2m',dei(i,k)
               !write(iulog,*) 'grid values of deff ice from optics file',g_d_eff
               !call endrun ('deff of ice exceeds limits')
            !endif
            ! for each cell interpolate to find weights and indices in g_d_eff grid.
            if (dlimited <= g_d_eff(1)) then
               k_d_eff = 2
               wd = 1._r8
               onemwd = 0._r8
            elseif (dlimited >= g_d_eff(n_g_d)) then
               k_d_eff = n_g_d
               wd = 0._r8
               onemwd = 1._r8 
            else
               do i_d_grid = 1, n_g_d
                  k_d_eff = i_d_grid
                  if(g_d_eff(i_d_grid) > dlimited) exit
               enddo
               wd = (g_d_eff(k_d_eff) - dlimited)/(g_d_eff(k_d_eff) - g_d_eff(k_d_eff-1))
               onemwd = 1._r8 - wd
            endif
            ! interpolate into grid and extract radiative properties
            do i_swband = 1, nswbands
               ext = wd*ext_sw_ice(k_d_eff-1,i_swband) + &
                 onemwd*ext_sw_ice(k_d_eff  ,i_swband) 
               ssa = wd*ssa_sw_ice(k_d_eff-1,i_swband) + &
                 onemwd*ssa_sw_ice(k_d_eff  ,i_swband) 
               asm = wd*asm_sw_ice(k_d_eff-1,i_swband) + &
                 onemwd*asm_sw_ice(k_d_eff  ,i_swband) 
               tau    (i_swband,i,k)=iciwpth(i,k) * ext
               tau_w  (i_swband,i,k)=iciwpth(i,k) * ext * ssa
               tau_w_g(i_swband,i,k)=iciwpth(i,k) * ext * ssa * asm
               tau_w_f(i_swband,i,k)=iciwpth(i,k) * ext * ssa * asm * asm
            enddo
         endif
      enddo
   enddo

   return
end subroutine get_ice_optics_sw   

!==============================================================================

subroutine get_liquid_optics_sw(state, pbuf, tau, tau_w, tau_w_g, tau_w_f)
   type(physics_state), intent(in) :: state
   type(pbuf_fld),      intent(in) :: pbuf(pbuf_size_max)

   real(r8),intent(out) :: tau    (nswbands,pcols,pver) ! extinction optical depth
   real(r8),intent(out) :: tau_w  (nswbands,pcols,pver) ! single scattering albedo * tau
   real(r8),intent(out) :: tau_w_g(nswbands,pcols,pver) ! asymetry parameter * tau * w
   real(r8),intent(out) :: tau_w_f(nswbands,pcols,pver) ! forward scattered fraction * tau * w

   real(r8), pointer, dimension(:,:) :: lamc, pgam, iclwpth
   real(r8), dimension(pcols,pver) :: kext
   integer i,k,swband,lchnk,ncol

   lchnk = state%lchnk
   ncol = state%ncol

   lamc => pbuf(i_lambda)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
   pgam => pbuf(i_mu)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
   iclwpth => pbuf(i_iclwp)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
   
   do k = 1,pver
      do i = 1,ncol
         if(lamc(i,k) > 0._r8) then ! This seems to be clue from microphysics of no cloud
            call gam_liquid_sw(iclwpth(i,k), lamc(i,k), pgam(i,k), &
                tau(1:nswbands,i,k), tau_w(1:nswbands,i,k), tau_w_g(1:nswbands,i,k), tau_w_f(1:nswbands,i,k))
         else
            tau(1:nswbands,i,k) = 0._r8
            tau_w(1:nswbands,i,k) = 0._r8
            tau_w_g(1:nswbands,i,k) = 0._r8
            tau_w_f(1:nswbands,i,k) = 0._r8
         endif
      enddo
   enddo

end subroutine get_liquid_optics_sw

!==============================================================================

subroutine liquid_cloud_get_rad_props_lw(state, pbuf, abs_od)
   type(physics_state), intent(in)  :: state
   type(pbuf_fld),      intent(in)  :: pbuf(pbuf_size_max)
   real(r8), intent(out) :: abs_od(nlwbands,pcols,pver)

   integer :: lchnk, ncol
   real(r8), pointer, dimension(:,:) :: lamc, pgam, iclwpth

   integer lwband, i, k

   abs_od = 0._r8

   lchnk = state%lchnk
   ncol = state%ncol

   lamc => pbuf(i_lambda)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
   pgam => pbuf(i_mu)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
   iclwpth => pbuf(i_iclwp)%fld_ptr(1,1:pcols,1:pver,lchnk,1)

   do k = 1,pver
      do i = 1,ncol
         if(lamc(i,k) > 0._r8) then ! This seems to be the clue for no cloud from microphysics formulation
            call gam_liquid_lw(iclwpth(i,k), lamc(i,k), pgam(i,k), abs_od(1:nlwbands,i,k))
         else
            abs_od(1:nlwbands,i,k) = 0._r8
         endif
      enddo
   enddo

end subroutine liquid_cloud_get_rad_props_lw
!==============================================================================

subroutine snow_cloud_get_rad_props_lw(state, pbuf, abs_od)
   type(physics_state), intent(in)  :: state
   type(pbuf_fld),      intent(in)  :: pbuf(pbuf_size_max)
   real(r8), intent(out) :: abs_od(nlwbands,pcols,pver)

   real(r8), pointer, dimension(:,:) :: dei, iciwpth
   real(r8) :: dlimited ! d limited to range dmin,dmax

   integer :: i,k,i_d_grid,k_d_eff,i_lwband,ncol,lchnk
   real :: wd, onemwd, absor

   abs_od = 0._r8

   ncol = state%ncol
   lchnk = state%lchnk

! note that this code makes the "ice path" point to the "snow path from CAM"
   iciwpth => pbuf(i_icswp)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
   dei     => pbuf(i_des)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
! note that this code makes the "ice path" point to the "snow path from CAM"

   do i = 1,ncol
      do k = 1,pver
         dlimited = dei(i,k) ! min(dmax,max(dei(i,k),dmin))
         ! if ice water path is too small, OD := 0
         if( iciwpth(i,k) < 1.e-80_r8 .or. dlimited .eq. 0._r8) then
            abs_od (:,i,k) = 0._r8
         !else if (dlimited < g_d_eff(1) .or. dlimited > g_d_eff(n_g_d)) then
         !   write(iulog,*) 'dlimited prognostic cldwat2m',dlimited
         !   write(iulog,*) 'grid values of deff ice from optics file',g_d_eff(1),' -> ',g_d_eff(n_g_d)
         !   !call endrun ('deff of ice exceeds limits')
         else
            ! for each cell interpolate to find weights and indices in g_d_eff grid.
            if (dlimited <= g_d_eff(1)) then
               k_d_eff = 2
               wd = 1._r8
               onemwd = 0._r8
            elseif (dlimited >= g_d_eff(n_g_d)) then
               k_d_eff = n_g_d
               wd = 0._r8
               onemwd = 1._r8 
            else
               do i_d_grid = 2, n_g_d
                  k_d_eff = i_d_grid
                  if(g_d_eff(i_d_grid) > dlimited) exit
               enddo
               wd = (g_d_eff(k_d_eff) - dlimited)/(g_d_eff(k_d_eff) - g_d_eff(k_d_eff-1))
               onemwd = 1._r8 - wd
            endif
            ! interpolate into grid and extract radiative properties
            do i_lwband = 1, nlwbands
               absor = wd*abs_lw_ice(k_d_eff-1,i_lwband) + &
                   onemwd*abs_lw_ice(k_d_eff  ,i_lwband)
               abs_od (i_lwband,i,k)=  iciwpth(i,k) * absor 
            enddo
         endif
      enddo
   enddo

end subroutine snow_cloud_get_rad_props_lw

!==============================================================================

subroutine ice_cloud_get_rad_props_lw(state, pbuf, abs_od)
   type(physics_state), intent(in)  :: state
   type(pbuf_fld),      intent(in)  :: pbuf(pbuf_size_max)
   real(r8), intent(out) :: abs_od(nlwbands,pcols,pver)

   real(r8), pointer, dimension(:,:) :: dei, iciwpth
   real(r8) :: dlimited ! d limited to range dmin,dmax

   integer :: i,k,i_d_grid,k_d_eff,i_lwband,ncol,lchnk
   real :: wd, onemwd, absor

   abs_od = 0._r8

   ncol = state%ncol
   lchnk = state%lchnk

   iciwpth => pbuf(i_iciwp)%fld_ptr(1,1:pcols,1:pver,lchnk,1)
   dei     => pbuf(i_dei)%fld_ptr(1,1:pcols,1:pver,lchnk,1)

   do i = 1,ncol
      do k = 1,pver
         dlimited = dei(i,k) ! min(dmax,max(dei(i,k),dmin))
         ! if ice water path is too small, OD := 0
         if( iciwpth(i,k) < 1.e-80_r8 .or. dlimited .eq. 0._r8) then
            abs_od (:,i,k) = 0._r8
         !else if (dlimited < g_d_eff(1) .or. dlimited > g_d_eff(n_g_d)) then
         !   write(iulog,*) 'dlimited prognostic cldwat2m',dlimited
         !   write(iulog,*) 'grid values of deff ice from optics file',g_d_eff(1),' -> ',g_d_eff(n_g_d)
         !   !call endrun ('deff of ice exceeds limits')
         else
            ! for each cell interpolate to find weights and indices in g_d_eff grid.
            if (dlimited <= g_d_eff(1)) then
               k_d_eff = 2
               wd = 1._r8
               onemwd = 0._r8
            elseif (dlimited >= g_d_eff(n_g_d)) then
               k_d_eff = n_g_d
               wd = 0._r8
               onemwd = 1._r8 
            else
               do i_d_grid = 2, n_g_d
                  k_d_eff = i_d_grid
                  if(g_d_eff(i_d_grid) > dlimited) exit
               enddo
               wd = (g_d_eff(k_d_eff) - dlimited)/(g_d_eff(k_d_eff) - g_d_eff(k_d_eff-1))
               onemwd = 1._r8 - wd
            endif
            ! interpolate into grid and extract radiative properties
            do i_lwband = 1, nlwbands
               absor = wd*abs_lw_ice(k_d_eff-1,i_lwband) + &
                   onemwd*abs_lw_ice(k_d_eff  ,i_lwband)
               abs_od (i_lwband,i,k)=  iciwpth(i,k) * absor 
            enddo
         endif
      enddo
   enddo

end subroutine ice_cloud_get_rad_props_lw

!==============================================================================

subroutine gam_liquid_lw(clwptn, lamc, pgam, abs_od)
  real(r8), intent(in) :: clwptn ! cloud water liquid path new (in cloud) (in g/m^2)?
  real(r8), intent(in) :: lamc   ! prognosed value of lambda for cloud
  real(r8), intent(in) :: pgam   ! prognosed value of mu for cloud
  real(r8), intent(out) :: abs_od(1:nlwbands)
  ! for interpolating into mu/lambda
  integer :: imu, kmu, wmu, onemwmu
  integer :: ilambda, klambda, wlambda, onemwlambda, lambdaplus, lambdaminus
  integer :: lwband ! sw band index
  real(r8) :: absc

  if (clwptn < 1.e-80_r8) then
    abs_od = 0._r8
    return
  endif

  if (pgam < g_mu(1) .or. pgam > g_mu(nmu)) then
    write(iulog,*)'pgam from prognostic cldwat2m',pgam
    write(iulog,*)'g_mu from file',g_mu
    call endrun ('pgam exceeds limits')
  endif
  do imu = 1, nmu
    kmu = imu
    if (g_mu(kmu) > pgam) exit
  enddo
  wmu = (g_mu(kmu) - pgam)/(g_mu(kmu) - g_mu(kmu-1))
  onemwmu = 1._r8 - wmu

  do ilambda = 1, nlambda
    klambda = ilambda
    if (wmu*g_lambda(kmu-1,ilambda) + onemwmu*g_lambda(kmu,ilambda) < lamc) exit
  enddo
  if (klambda < 1 .or. klambda > nlambda) call endrun('lamc  exceeds limits')
  lambdaplus = wmu*g_lambda(kmu-1,klambda  ) + onemwmu*g_lambda(kmu,klambda  )
  lambdaminus= wmu*g_lambda(kmu-1,klambda-1) + onemwmu*g_lambda(kmu,klambda-1)
  wlambda = (lambdaplus - lamc) / (lambdaplus - lambdaminus)
  onemwlambda = 1._r8 - wlambda

  do lwband = 1, nlwbands
     absc=     wlambda*    wmu*abs_lw_liq(kmu-1,klambda-1,lwband) + &
           onemwlambda*    wmu*abs_lw_liq(kmu-1,klambda  ,lwband) + &
               wlambda*onemwmu*abs_lw_liq(kmu  ,klambda-1,lwband) + &
           onemwlambda*onemwmu*abs_lw_liq(kmu  ,klambda  ,lwband)

     abs_od(lwband) = clwptn * absc
  enddo

  return
end subroutine gam_liquid_lw

!==============================================================================

subroutine gam_liquid_sw(clwptn, lamc, pgam, tau, tau_w, tau_w_g, tau_w_f)
  real(r8), intent(in) :: clwptn ! cloud water liquid path new (in cloud) (in g/m^2)?
  real(r8), intent(in) :: lamc   ! prognosed value of lambda for cloud
  real(r8), intent(in) :: pgam   ! prognosed value of mu for cloud
  real(r8), intent(out) :: tau(1:nswbands), tau_w(1:nswbands), tau_w_f(1:nswbands), tau_w_g(1:nswbands)
  ! for interpolating into mu/lambda
  integer :: imu, kmu, wmu, onemwmu
  integer :: ilambda, klambda, wlambda, onemwlambda, lambdaplus, lambdaminus
  integer :: swband ! sw band index
  real(r8) :: ext, ssa, asm

  if (clwptn < 1.e-80_r8) then
    tau = 0._r8
    tau_w = 0._r8
    tau_w_g = 0._r8
    tau_w_f = 0._r8
    return
  endif

  if (pgam < g_mu(1) .or. pgam > g_mu(nmu)) then
    write(iulog,*)'pgam from prognostic cldwat2m',pgam
    write(iulog,*)'g_mu from file',g_mu
    call endrun ('pgam exceeds limits')
  endif
  do imu = 1, nmu
    kmu = imu
    if (g_mu(kmu) > pgam) exit
  enddo
  wmu = (g_mu(kmu) - pgam)/(g_mu(kmu) - g_mu(kmu-1))
  onemwmu = 1._r8 - wmu

  do ilambda = 1, nlambda
     klambda = ilambda
     if (wmu*g_lambda(kmu-1,ilambda) + onemwmu*g_lambda(kmu,ilambda) < lamc) exit
  enddo
  if (klambda < 1 .or. klambda > nlambda) call endrun('lamc  exceeds limits')
     lambdaplus = wmu*g_lambda(kmu-1,klambda  ) + onemwmu*g_lambda(kmu,klambda  )
     lambdaminus= wmu*g_lambda(kmu-1,klambda-1) + onemwmu*g_lambda(kmu,klambda-1)
     wlambda = (lambdaplus - lamc) / (lambdaplus - lambdaminus)
  onemwlambda = 1._r8 - wlambda

  do swband = 1, nswbands
     ext =     wlambda*    wmu*ext_sw_liq(kmu-1,klambda-1,swband) + &
           onemwlambda*    wmu*ext_sw_liq(kmu-1,klambda  ,swband) + &
               wlambda*onemwmu*ext_sw_liq(kmu  ,klambda-1,swband) + &
           onemwlambda*onemwmu*ext_sw_liq(kmu  ,klambda  ,swband)
     ! probably should interpolate ext*ssa
     ssa =     wlambda*    wmu*ssa_sw_liq(kmu-1,klambda-1,swband) + &
           onemwlambda*    wmu*ssa_sw_liq(kmu-1,klambda  ,swband) + &
               wlambda*onemwmu*ssa_sw_liq(kmu  ,klambda-1,swband) + &
           onemwlambda*onemwmu*ssa_sw_liq(kmu  ,klambda  ,swband)
     ! probably should interpolate ext*ssa*asm
     asm =     wlambda*    wmu*asm_sw_liq(kmu-1,klambda-1,swband) + &
           onemwlambda*    wmu*asm_sw_liq(kmu-1,klambda  ,swband) + &
               wlambda*onemwmu*asm_sw_liq(kmu  ,klambda-1,swband) + &
           onemwlambda*onemwmu*asm_sw_liq(kmu  ,klambda  ,swband)
     ! compute radiative properties
     tau(swband) = clwptn * ext
     tau_w(swband) = clwptn * ext * ssa
     tau_w_g(swband) = clwptn * ext * ssa * asm
     tau_w_f(swband) = clwptn * ext * ssa * asm * asm
  enddo

  return
end subroutine gam_liquid_sw

!==============================================================================

end module cloud_rad_props
